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REDUCTION OF LARGE DYNAMICAL SYSTEMS BY MINIMIZATION OF 

EVOLUTION RATE 

SHARATH S. GIRIMAJI* 

Abstract. Reduction of a large system of equations to a lower- dimensional system of similar dynamics 
is investigated. For dynamical systems with disparate timescales, a criterion for determining redundant 
dimensions and a general reduction method based on the minimization of evolution rate are proposed. 
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1. Introduction. The macroscopic behavior of many complex systems (with a large number of degrees 
of freedom or scales) is largely insensitive to the details of the microscopic features. The macroscopic behavior 
can then be nearly exactly described by a simpler system with very few degrees of freedom or scales. In 
this article, a general method for reducing the description of the macroscopic behavior of large systems 
is presented. The two specific examples considered in this paper are of importance in the modeling and 
computation of turbulent combustion: the reduction of complex chemical kinetics and algebraic modeling of 
Reynolds stresses. 

Consider an autonomous dissipative dynamical system with one attracting fixed point: 

(1.1) z = g(z), where z = (zi, z 2 , • • • z n ). 

It is assumed that z is suitably scaled and nondimesionalized. If the system has disparate timescales, the 
solution exhibits a typical three-stage behavior, (i) Initial- condition dependent initial transient stage which 
lasts until all the small-timescale fast processes are exhausted; (ii) the intermediate slow-manifold stage in 
which the solutions ‘bunch’ together in a lower- dimensional phase space as the slow processes dominate; and 
(iii) the final equilibrium state. In the slow- manifold stage, the degrees of freedom of the system can be 
reduced if the relationship between fast and slow variables can be found. In a nonlinear dynamical system 
it is difficult to characterize the slow manifold accurately and the current practice is to locally linearize the 
equations. We will use the example of a general two- variable linear system with one large (“1) and one 
small (— e) eigenvalue in our analysis: 

x = — x(cos 2 6 + e sin 2 0) + y sin 0 cos 0(1 — e), 

(1.2) y = x sin 6 cos 0(1 — e) — y(sin 2 0 + 5 cos 2 9). 

The directions corresponding to the two eigenvalues are 

(1.3) x = —y cot 0 for A — — 1; x = +ytan0 for A = — e. 

In a linear system, the slow manifold lies in a space spanned by the eigenvectors corresponding to small 
(in magnitude) eigenvalues. Maas and Pope [1] determine the slow- manifold by requiring it to be orthogonal 
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to the principal directions corresponding to the large (negative) eigenvalues of the Jacobian (MP model). In 
the example considered the slow-manifold is 


(1.4) 


x(slow- manifold) = f y tan 0. 


The determination of the eigenvalues and eigen-directions of the Jacobian and subsequently the slow manifold 
can be prohibitively expensive in a practical nonlinear problem involving several dozen variables. The popular 
and inexpensive steady-state approximation (SSA) method of reduction involves setting the rate of change 
of fast variables to zero to obtain algebraic relations between fast and slow variables. In the sample problem, 

x (slow-manifold) = +y tan 9 (l - e sec 2 9 + e 2 sec 2 9 tan 2 9 + 0(e 3 )) , if x is fast variable, 

(1.5) y(slow-manifold) = +xcot0 (l - e esc 2 9 + e 2 esc 2 9 cot 2 9 + 0(e 3 )) , if y is fast variable. 

Judicious choice of the steady state variable is crucial and this is usually made from a priori knowledge. 
Even with a proper choice, the error incurred is first order in the eigenvalue ratio. For further details about 
these methods and other important developments the reader is referred to [2], [3] and [4], 

Our objectives are to (i) devise a general criterion for selecting fast variables and (ii) develop a compu- 
tationally viable reduction procedure of higher (than one) order accuracy in eigenvalue ratio. 

2. Characterization of Slow Variables. Fast variables are those whose evolution is dominated by 
the large eigenvalues of the Jacobian and the slow variables by the small eigenvalues. In fact, close to the 
slow-manifold, fast variable evolution can be slower than that of slow variables rendering evolution rate an 
unsuitable selection criterion. A better gauge of the timescale of a variable (z, ) is the convergence rate of its 
value between two neighboring trajectories. The difference in z t values between two neighboring trajectories, 
Szi, evolves according to 



Based on this, two approximate measures of the dissipativeness of individual variables are proposed for the 
selection criterion: 


( 2 . 2 ) 



or Di = 



In general, large magnitude of Di implies fast variable and slow variables are characterized by small magni- 
tudes. When the physical and eigen variables coincide, it is easy to see that D t are indeed the eigenvalues. 
While the second estimate is likely to be more accurate in general, the first measure may be adequate when 
the off-diagonal terms in the Jacobian are relatively small. 


3. Reduction of the Dynamical System. Each state in phase-space is associated with a residence 
timescale which is proportional to the amount of time a solution trajectory resides in an infinitesimal neigh- 
borhood of that point and inversely proportional to the local evolution rate (V = vCCiU HifU)- Solutions 
tend to stagnate near long-timescale states and pass quickly through short-timescale states. The probability 
of finding a solution, of unknown time lapse and arbitrary initial condition, at a given state along its trajec- 
tory is proportional to residence timescale of that state. This leads to the premise of our reduction procedure. 
Given the values of slow variables, an arbitrary solution trajectory is most likely to be found at the state with 
the longest residence time, i.e., smallest evolution rate. This proposal is as valid for non-linear systems as 
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it is for linear systems and bears semblance to the maximum likelihood estimator of mathematical statistics 
and the non- equilibrium potential concept of [5]. When not conditioned by any slow- variable values, this 
criterion selects the equilibrium state as the most likely state of the system. If we represent the retained 
(slow) variables in z by the vector y and the discarded (fast) variables by the vector x, the statement of our 
proposal can be written aS; 

n 

(3.1) x(y) « minF(x : y), where F = V 2 = ^ ff,g t . 

x »=i 


In the present example 


(3.2) 


F(x, y) 


dx 

2 

dy 

dt 

+ 

dt 


= a: 2 (cos 2 6 + e 2 sin 2 6) + y 2 ( sin 2 6 + e 2 cos 2 0) — 2xys'm6cos0(l — £ 2 ). 


The value of x that minimizes the evolution rate is 


(3.3) 


sin 9 cos 0(1 — e 2 ) 

^ cos 2 9 + e 2 sin 2 9 

= ytan0[l — (1 + tan 2 0)e 2 + 0(e 4 )]. 


When the eigen and physical variables coincide, there is no error involved; otherwise, it is of order e 2 which 
makes this method more accurate than the SSA method. The fast variable selection criterion given above 
will ensure that 0 < n/4. 

For even better accuracy, minimization of the evolution rates of higher-order derivatives of the variables 
is proposed. Differentiation has the effect of separating the timescales so that the timescales of x and y are 
farther apart than the those x and y. The higher the order of the derivative, the greater is the separation 
between the slow and fast variable (derivatives). In matrix algebra, this fact is the basis of the power method 
of separating small and large eigenvalues. Our proposal for higher-order accuracy is 

n d m z d™ z 

(3.4) x (y) « minf""(* : y); where F"*( z) = £ 

i=\ 

In our sample problem, 


F m (x, y) = a; 2 (cos 2 6 + e 2rn sin 2 9) + y 2 (sin 2 9 + e 2m cos 2 9) 

(3.5) — 2xy sin0cos0(l — e 2m ), 

which leads to 

(3.6) x(y) = y tan0[l - (1 + tan 2 0)e 2m + 0(e 4m )]. 

In practical problems, minimization of F 1 may be adequate as the error involved in local linearization 
may be of order e 2 . More than one minima of the evolution rate may be encountered in the proximity of 
other (unstable) fixed points. (This will also be the case with SSA method.) Then, the minimum that has ah 
negative eigenvalues (attractor) must be chosen. The computational effort is still much smaller than the MP 
method, where eigenvalues and eigenvectors are evaluated at all points in the fast phase-space as opposed 
to eigenvalues only at a handful of minima in the present method. 
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4. Chemical Kinetics Reduction. The phenomenon of turbulent combustion encompasses a wide 
spectrum of length and time scales. Typically, the thermo- chemical timescales representing chemical reactions 
span a wider range than those of turbulent advection and molecular diffusion. In the interest of computational 
efficiency, it is desirable to eliminate chemical reactions and species whose characteristic timescales are smaller 
than the fluid timescales. 


Perfectly-stirred reactor 



Fig. 4.1. State-space evolution of H and OH as a function of temperature . Different dashed and dotted lines correspond 
to results from full kinetics simulations from different initial conditions. The solid line represents the slow-manifold model 
prediction with temperature as the only retained variable. 


Consider chemical reactions among IV s chemical species involving Af r reactions in an adiabatic, isobaric, 
well-stirred reactor. The species (mass fraction 1^) evolve according to 


(4.1) 


BY Nr 


J = 1 


Pij 

P<*3 


AjT bj 


exp 


(-&) 8 (&) 


where: W t is the molecular weight of species t; fa is the stoichiometric coefficient of species i in the j-th 
reaction; aj = 5Z»=i Pij > 4?> all( l Ej are the Arrehenius constants; and, 7Z is the universal gas constant 
Temperature and density are obtained from 


(4.2) 



Ns 

yi A/ii, o 

i= 1 


dYj 
~dt ’ 


p(p,T,Yi) 


P 

nr 
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We consider the evolution of a H 2 /0 2 mixture, in which the atomic mass fractions of hydrogen and oxygen 
are 0.15 and 0.85 respectively, from various initial conditions. The constant total enthalpy is 1.104 x 10 7 
KJ and the constant pressure is 2 atm. The detailed 14-step chemical mechanism involves temperature (T) 
and six species (0 2 , O, H 2 , H, OH, and H 2 0). The equilibrium state is given by: T = 3506 K : Y]] = 0.024: 
5"h 2 o = 0.556; y OH = 0.201; Y 0 = 0.086; y Ha = 0.060 and Yq 2 = 0.070. 


Perfectly— stirred reactor 



1 2 3 4. 

Temperature 

Fig. 4.2. State-space evolution of H and OH as a function of temperature. Solid line - case IC5 with full kinetics ; dotted 
line - reduced model with temperature as the only retained variable; dot-dash line - reduced model with temperature and H 2 0 
as the retained variables. 


Evolution of the species undergoing reactions according to the detailed chemical kinetics from six arbi- 
trary initial conditions is calculated numerically using a second-order Runge-Kutta procedure. The evolution 
trajectories of H and OH are shown in Figure 4.1. A one-dimensional slow manifold model is constructed 
with temperature as the only retained variable: 

(4-3) Y (T) % min 

Y 

The mass fractions of H and OH predicted by the one- dimensional model is also shown in Figure 4.1. The 
agreement is quite adequate. Next we construct a two dimensional manifold with water mass fraction as 
a retained variable along with temperature. Comparison between the two-dimensional model and detailed 
kinetics is made for the case IC5 in Figure 4.2. The one-dimensional model values are also shown for 
comparison. The agreement between the two-dimensional model and detailed kinetics is, in general, excellent. 
These results demonstrate that the present method is an accurate and computationally viable option for 


N t 


E n n dTdT 


i=l 
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reducing chemical kinetics. 


5. Algebraic Reynolds Stress Modeling. In the field of turbulence modeling, closure of Reynolds 
stress transport equations is the lowest level of sophistication at which models can be developed in a sys- 
tematic manner from the governing Navier-Stokes equations. This entails solving seven modeled transport 
equations for turbulence variables five for Reynolds stress anisotropy components b{j and one each for 
turbulent kinetic energy K and dissipation e 

(5.1) ^ - L\b mn S mn ) 

2 

-\-L2Sij + L^^bikSjk *+■ bjkSik ~~ ^bim^lm^ij ) 

+L^(bikW jk + bjkWik) 

dK 

- = - KbijSij -e 
§ = -C el eb l3 S l3 - Cj-. 

The Einstein summation convention is used. In the above equations L’s and C’s are model coefficients, and 
Sij and Wij are the mean flow strain and rotation rates. For the sake of brevity, the coefficient values used 
are not given here and the reader is referred to [6]. 
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Fig. 5.1. Homogeneous turbulence. Evolution 0/612 and 611 os a function of u>. Solid line - reduced model with K and e 
as the retained variables ; dashed and dotted lines - full Reynolds stress model calculations from various initial conditions. 
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Researchers have long sought to reduce the seven-equation model to a computationally affordable two- 
equation (for K and e) model. The Reynolds stress anisotropies are the discarded variables for which we 
seek algebraic expressions: hence the name, algebraic stress model. The evolution rate of this system in 
state- space is given by, 

/c 2 ) F = V 2 — | de 

' ' dt dt dt dt dt dt 

The algebraic stress model assumption for non-equilibrium turbulent flows is 

(5.3) b (AT, e) = min F(h : K, e). 

b 

For two-dimensional mean flows, the minimization is performed in a three-dimensional discarded variable 
state space of &n, 612 and 633 . 

We test the model in homogeneous shear turbulence in which the only non-zero values of mean strain 
and rotation tensor components are 

(5.4) S\ 2 = 521 = R 12 = — W 21 = — , 

where 5 is a constant. The seven-equation model calculations for various 5 values are shown in Figure 5.1. 
Each anisotropy is plotted as a function of the self-similarity (retained) variable e/SK. The equilibrium 
values are: u> = 0.166; bn = 0.204; 612 = —0.157; 622 = —0.149; and 633 = —0.055. The various solution 
trajectories exhibit a well-defined slow manifold behavior. The algebraic model, also shown in the figure, 
does an excellent job of reproducing the behavior of the full set. This algebraic model has been tested in a 
variety of other homogeneous flows with good success [7]. 

6 . Conclusion. In summary, we have developed a general procedure for reducing dimensionality of 
autonomous systems with disparate timescales. The premise of the proposal is that the solution trajectories 
from unknown initial conditions are most likely to be found near the ‘bottle- necks’ in state-space which are 
characterized by small evolution rates. The location of these bottle-necks on the path to equilibrium can be 
found by simply minimizing the evolution velocity subject to constraints (given values of the slow or retained 
variables) . 
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